home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / lang / c-part1 / 4052 < prev    next >
Encoding:
Text File  |  1996-08-05  |  6.3 KB  |  210 lines
  1. Path: news.vanderbilt.edu!news
  2. From: "Marcus H. Mendenhall" <mendenm@ctrvax.vanderbilt.edu>
  3. Newsgroups: comp.lang.c
  4. Subject: Re: help with pi algorithm
  5. Date: Thu, 01 Feb 1996 10:17:56 +0000
  6. Organization: Vanderbilt University, Nashville, TN, USA
  7. Message-ID: <31109354.13A1@ctrvax.vanderbilt.edu>
  8. References: <0fc_9601240111@csource.blaze.net.au> <4ehfmb$1kn@guava.epix.net>
  9. NNTP-Posting-Host: 129.59.235.1
  10. Mime-Version: 1.0
  11. Content-Type: text/plain; charset=us-ascii
  12. Content-Transfer-Encoding: 7bit
  13. X-Mailer: Mozilla 2.0b6a (Macintosh; I; PPC)
  14.  
  15. Here is a piece of code I wrote which portably computes PI.  I have computer 1,000,000 
  16. digits with a PowerMacintosh (took about a day, if I remember correctly).  Notice that the 
  17. iostreams stuff can be easily removed and replaced with standard C-style i/o if you are 
  18. using c and not c++.  To compute lots of digits on a PC running in 16-bit mode, you will 
  19. have to declare some things huge.  Good Luck.
  20.  
  21.  
  22. #include <iostream>
  23. #include        <iomanip>
  24. #include        <stdlib.h>
  25. #include        <fp.h>
  26. #include        <profiler.h>
  27. #include        <time.h>
  28.  
  29. unsigned long longmult(long dwords, unsigned long mult, unsigned long data[]);
  30. int longdiv(long dwords, long divisor, unsigned long data[], int start);
  31. int longadd(long dwords, unsigned long src[], unsigned long dst[]);
  32. int longsub(long dwords, unsigned long src[], unsigned long dst[]);
  33.  
  34. void main(void)
  35. {
  36.         unsigned long *sum, *p8, *p57, *p239, *scr;
  37.         clock_t start, stop;
  38.         
  39.         long lwords;
  40.         int zt8=0, zt57=0, zt239=0, index;
  41.         
  42.         ofstream pidata("pi_digits.txt");
  43.         
  44.         cout << "Enter number of longwords: ";
  45.         cin >> lwords;
  46.         
  47.         scr = new unsigned long[lwords];
  48.         sum = new unsigned long[lwords];
  49.         p8 = new unsigned long[lwords];
  50.         p57 = new unsigned long[lwords];
  51.         p239 = new unsigned long[lwords];
  52.         
  53.         if(!p239) {
  54.                 cout << "Not enough memory..." << endl;
  55.                 exit(1);
  56.         }
  57.         
  58.         start=clock();
  59.         
  60. #if     __profile__
  61.         OSErr err;
  62.         
  63.         err=ProfilerInit(collectDetailed, bestTimeBase, 100, 100);
  64.         if(err) {
  65.                 cout << "Profiler open error " << err << endl;
  66.                 exit(1);
  67.         }
  68. #endif
  69.  
  70.         for(int i=0; i<lwords; i++) {
  71.                 sum[i]=0; p8[i]=0; p57[i]=0; p239[i]=0;
  72.         }
  73.         
  74.         p8[0]=24; p57[0]=8; p239[0]=4;
  75.         longdiv(lwords, 8, p8, 0);
  76.         longdiv(lwords, 57, p57, 0);
  77.         longdiv(lwords, 239, p239, 0);
  78.         
  79.         index=1;
  80.         while(zt8 < lwords) {
  81.                 for(int j=0; j<lwords; j++) scr[j]=p8[j];
  82.                 if (zt57 < lwords) longadd(lwords, p57, scr);
  83.                 if (zt239 < lwords) longadd(lwords, p239, scr);
  84.                 longdiv(lwords, index, scr, zt8 > 0 ? zt8-1 : 0);
  85.                 if (index & 2) longsub(lwords, scr, sum);
  86.                 else longadd(lwords, scr, sum);
  87.                 zt8=longdiv(lwords, 8*8, p8, zt8);
  88.                 if(zt57 < lwords) zt57=longdiv(lwords, 57*57, p57, zt57);
  89.                 if(zt239 < lwords) zt239=longdiv(lwords, 239*239, p239, zt239);
  90.                 index += 2;
  91.         }
  92.  
  93.         stop = clock();
  94.         
  95. #if     __profile__
  96.         err=ProfilerDump("\pPi_Profile.dat");
  97.         if(err) {
  98.                 cout << "Profiler dump error " << err << endl;
  99.         }
  100.         ProfilerTerm();
  101. #endif
  102.  
  103.         cout << "Index = " << index << endl;
  104.         
  105.         for(int i=0; i<lwords; i++) {
  106.                 pidata << setw(9) << right << setfill('0') << sum[i] << " ";
  107.         }
  108.         pidata << endl;
  109.         pidata << "Time = " << (float)(stop-start)/CLOCKS_PER_SEC << " seconds" << endl;
  110. }
  111.  
  112.  
  113. const static double shift=1000000000.0, shifti=1.0/(1000000000.0),
  114.         epsi = 0.5/(1000000000.0);
  115.  
  116. /* longdiv divides an arbitrary-precision number represented with long ints by divisor.
  117.         there are 31 sig bits / long int (to leave room for overflows)
  118. */
  119. int longdiv(long dwords, long divisor, unsigned long data[], int start)
  120. {
  121.         register double d, dinv, eps2, carry;
  122.         register double d1, fd, rshift=shift;
  123.         register unsigned long r1;
  124.         
  125.         int zt=0, i;
  126.         
  127.         d=divisor;
  128.         dinv=1.0/d;
  129.         eps2=0.5*dinv;
  130.         
  131.         carry=0.0;
  132.         
  133.         for(i=start; i<dwords; i++) if(data[i]) break; // find non-zero start
  134.         
  135.         zt=i;
  136.         
  137.         for(; i<dwords; i++) {
  138. #ifndef GENERATINGPOWERPC
  139.                 d1=data[i]+carry*rshift;
  140.                 r1=(int)(d1*dinv+eps2);
  141.                 carry=d1-r1*d;
  142.                 data[i]=r1;
  143. #else
  144.                 fd=data[i];
  145.                 d1=fd*rshift;
  146.                 r1=(int)((fd+carry)*dinv+eps2);
  147.                 carry=(carry-r1*d)*rshift+d1;
  148.                 data[i]=r1;
  149. #endif
  150.  
  151.         }
  152.         return zt;
  153. }
  154.  
  155. /* longmult multiplies an arbitrary-precision number represented with long ints by mult.
  156.         there are 31 sig bits / long int (to leave room for overflows)
  157.         return value is carry out of highest word.
  158. */
  159. unsigned long longmult(long dwords, unsigned long mult, unsigned long data[])
  160. {
  161.         double carry;
  162.         double m;
  163.         
  164.         carry=0;
  165.         m=mult;
  166.         
  167.         for(int i=dwords-1; i>=0; i--) {
  168.                 double d1;
  169.                 d1=(data[i]*m)+carry;
  170.                 carry=(unsigned long)(d1*shifti+epsi); /* I hope converting to int 
  171. truncates */
  172.                 data[i]=(int)(d1-carry*shift);
  173.         }       
  174.         return carry;
  175. }
  176.  
  177. int longadd(long dwords, unsigned long src[], unsigned long dst[])
  178. {
  179.         int carry=0;
  180.         unsigned int sum;
  181.         for(int i=dwords-1; i>=0; i--) {
  182.                 sum=src[i]+dst[i]+carry;
  183.                 if(sum > 1000000000) {
  184.                         dst[i]=sum -1000000000;
  185.                         carry=1;
  186.                 } else {
  187.                         dst[i]=sum;
  188.                         carry=0;
  189.                 };
  190.         }
  191.         return carry;
  192. }
  193.  
  194. int longsub(long dwords, unsigned long src[], unsigned long dst[])
  195. {
  196.         int carry=0;
  197.         int sum;
  198.         for(int i=dwords-1; i>=0; i--) {
  199.                 sum=(int)dst[i]-(int)src[i]-carry;
  200.                 if(sum < 0) {
  201.                         dst[i]=sum + 1000000000;
  202.                         carry=1;
  203.                 } else {
  204.                         dst[i]=sum;
  205.                         carry=0;
  206.                 }
  207.         }
  208.         return carry;
  209. }
  210.